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highlights  in  the  following. 

Multi-reference  coupled-cluster  theory. 

One  objective  was  to  develop  our  proposed  multi-reference  coupled-cluster  method,  MIP/MEA-EOM-CCSD,  which  means  multi¬ 
ionization  potential,  multi-electron  attached,  equation  of  motion  coupled-cluster  (CC)  theory.  For  the  double  IP  and  double  EA 
(DIP/DEA)  case,  the  basic  idea  is  that  by  doing  a  CC  calculation  for  a  doubly  ionized  or  doubly  electron  attached  system,  the 
full  multi-reference  character  for  the  most  common,  2x2  problem  can  be  easily  introduced.  Then,  one  will  have  4  orbitals  that 
can  be  occupied  by  2  electrons  in  any  possible  way,  as  shown  below. 

Insert  fig  1 

So  for  cases  where  the  single  reference  (SR-CC)  is  inadequate,  because  CO  is  far  from  90%,  the  fact  that  we  can  obtain  the 
mixing  values,  C2,  Cl,  and  C3  from  the  EOM-CC  equations  enables  us  to  apply  this  method  easily  and  in  a  straight-forward 
way.  The  method  has  several  attractive  features  that  make  the  approach  worth  pursuing.  These  include  the  following: 

The  ansatz  has  a  global  extensive  part,  exp<^(TA(n±2) )  |n±2>  based  on  a  closed  shell  n  ±  2  reference  vacuum  and  a  local 
correlation  intensive  part,  whose  wave  function  is  |  |<^>=RA(n±2)  exp^(TA(n±2) )  |n±2>,  where  RA(n±2)  is  the  Cl-like  right-hand 
eigenvector  in  EOM-CC.  This  has  the  advantage  that  instead  of  asking  a  fully  extensive  MR-CC  method  to  account  for  dynamic 
and  non-dynamic  correlation  effects,  the  intensive  part  allows  one  to  target  the  usually  local  multi-reference  behavior. 

As  long  as  the  |n±2>  reference  is  a  closed-shell,  the  target  states  are  automatically  spin-eigenfunctions.  The  guarantee  of 
a  spin-eigenfunction  is  not  always  achievable  in  most  SR-CC  or  MR-CC  methods. 

The  DIP/DEA-EOM-CC  wavefunction  is  operationally  single  reference  making  it  as  easy  to  apply  as  single-reference  CC, 
with  no  decisions  for  the  user  but  basis  set,  level  of  correlation,  and  a  choice  of  the  one  spatial  active  orbital  to  doubly  occupy  in 
the  n  +  2  vacuum  (usually  the  LUMO),  or  to  un-occupy  in  the  n  -  2  vacuum  (usually  the  HOMO). 

The  DIP/DEA-EOM-CC  is  invariant  to  active  orbital  rotations  by  virtue  of  the  active  orbitals  being  in  either  the  occupied  or 
the  unoccupied  space. 

Multiple  states  can  be  obtained  from  the  EOM  matrix  diagonalization  providing  excited  states  as  well  as  the  ground  state. 

So  in  the  event  that  other  occupied  orbitals  in  DIP  or  virtuals  in  DEA  interact  strongly  with  the  chosen  active  orbital,  then  that 
solution  occurs  as  well,  and  will  appear  as  one  of  the  eigenvectors.  This  helps  to  confirm  the  particular  orbitals  that  manifest 
MR  character  and  if  there  are  more  than  two,  might  suggest  a  subsequent  three  (TIP/TEA)  or  four  (QIP/QEA)  calculation.  The 
price  paid  for  these  attractive  features  is  that  each  stage  of  the  calculation  has  to  be  converged:  first  the  SR-CC  solution  for  the 
n±2  closed  shell  system,  then  the  DEA/DIP-EOM  solution  itself.  Because  of  orbital  dependence  in  these  calculations,  the  DIP 
solution  in  particular  can  sometimes  be  difficult  to  converge. 

An  accurate  study  of  the  singlet-triplet  separations  in  di-radical  molecules  depends  upon  a  correct  description  of  the  multi¬ 
reference  character.  The  problem  is  the  prototype  for  2  electrons  in  two  orbitals,  I  and  A,  where  the  four  determinants  lallJ, 
AaAR>,  laAIJ  and  IIJAa,  could  have  major  weight  in  a  description  of  the  singlet  state,  that  is  not  recovered  by  a  single  reference 
calculation  based  upon  lallJ.  The  triplet  can  often  be  adequately  described  by  a  single  reference  description  based  upon  laAa, 
which  is  equivalent  to  the  triplet  coupling  of  the  two  laAli  determinants.  But  the  singlet  has  no  alternative  solution. 

Consequently,  we  applied  our  DIP/DEA  MR  methods  for  this  two  electron  two-orbital  problem. 

The  prototype  for  all  such  diradicals  is  methylene,  while  more  complicated  systems  like  the  three  forms  of  benzyne,  tri¬ 
methylene  methane  (TMM),  and  the  interesting  molecule  (CO)4  were  also  studied  in  this  paper  dedicated  to  Isaiah  Shavitt  [1]. 
See  below  for  an  illustration  of  the  benzyne  isomers. 

Ortho-,  meta-  and  para-benzyne  singlet-triplet  splittings  (in  kcal/mol)1. 


cc-pVDZ  cc- 

■pVTZ 

cc-pVQZ  3-4 

Extrap 

Ortho 

CCSD  28.6  30 

.6  31 

.3  31.8 

CCSD(T)  33.8 

37.0 

37.4  38.0 

^CCSD(T)  33.4 

36.6 

36.9  37.1 

DEA-EOM-CCSD 

33.5 

36.1  36.9 

37.5 

DEA-STEOM-CCSD  30.6  34.3  34.2  34.0 


DIP-EOM-CCSD 

37.5 

42.6 

46.3  49.0 

Experiment-ZPE-core 

37.4±0.3 

37.5±0.5 

Meta 

CCSD  9.9  10.5 

10.7 

10.8 

CCSD(T)  20.4 

22.1 

22.4 

22.7 

^CCSD(T)  20.9 

23.1 

20.9 

19.3 

DEA-EOM-CCSD 

17.1 

18.3 

18.4  18.5 

DEA-STEOM-CCSD  14.0  15.3  18.4  20.6 

DIP-EOM-CCSD 

18.0 

19.7 

18.3  17.3 

Experiment-ZPE-core 


20.5±0.3 


Para 

CCSD  -17.8  -19.2  -19.5  -19.7 

CCSD(T)  4.0  3.6  3.5  3.4 
^CCSD(T)  1.9  2.3  2.4  2.5 
DEA-EOM-CCSD  3.0  3.4  3.4  3.4 
DEA-STEOM-CCSD  4.9  4.6  4.3  4.0 
DIP-EOM-CCSD  3.9  4.4  N/C 
Experiment-ZPE-core  4.4±0.3 

2.7±0.4 

Our  answers  are  good,  but  there  is  still  an  issue  with  the  optimum  choice  of  the  orbitals  to  use  in  such  calculations.  In  some 
cases,  like  in  methylene,  the  DIP-EOM-CCSD  method  tends  to  diverge  when  basis  sets  are  very  large.  In  other  words,  the 
underlying  description  of  the  di-anion  suffers  from  an  admixture  of  the  continuum.  The  same  issue  pertains  to  the  lack  of 
convergence  in  para-benzyne.  This  does  not  happen  for  the  DEA-EOM  based  upon  the  double  cation.  Even  for  the  DIP-EOM 
this  divergence  does  not  happen  for  the  di-radical  isomers  TTM  or  (CO)4,  but  this  is  indicative  of  a  problem  that  needs  to  be 
solved  to  develop  a  truly  robust  method.  Orbital  optimization  would  be  one  approach;  constraining  potentials  another. 

Multiplets  Transition  Metal  and  Template  Orbitals 

Transition  metal  multiplets  are  a  prototype  of  another  type  of  multi-reference  problem.  The  five  d-orbital  degeneracy  plus  the 
close  lying  4s  orbital  leads  to  many  different  states  (multiplets)  depending  upon  their  electron  occupancy.  Consequently,  the 
small  energy  separations  among  these  mutiplets  are  a  sensitive  measure  of  the  ability  of  a  quantum  chemical  method  to  handle 
complicated  degeneracies  and  spin  states.  We  addressed  this  issue  and  applied  high-level  single  reference  coupled-cluster 
methods  to  define  a  benchmark.  At  the  CCSDT-3  level  we  obtained  excellent  results  [2]  and  demonstrated  that  the  results  were 
insensitive  to  a  variety  of  different  orbital  choices,  including  Brueckner  orbitals  and  fractionally  occupied  SCF  orbitals.  The  latter 
in  particular  is  our  new,  ‘template  orbital’  approach  that  offers  many  attractive  aspects. 

By  fractionally  occupying  the  orbitals  in  an  SCF  calculation,  like  describing  a  state  with  up  to  6  d  electrons  by  putting  3/5  of  an 
electron  in  each  spatial  orbital,  and  then  doing  an  SCF  calculation  to  define  the  set  of  orbitals  to  use  in  subsequent  CC 
calculations.  This  set  of  ‘template’  orbitals  is  then  occupied  to  accommodate  any  occupation  of  d  and  s  electrons  that  can  occur 
in  the  atom’s  multiplets.  Thus  one  set  of  orbitals  is  used  for  all  multiplet  states  that  correspond  to  different  occupancies  for  the 
3d4s  orbitals.  Such  orbitals  are  not  variationally  optimum  as  they  would  be  if  determined  from  separate  SCF.  calculations,  but 
unlike  the  latter,  which  would  constitute  a  different  calculation  for  every  multiplet  state,;  our  template  approach  exploits  the  fact 
the  CCSD  and  beyond  will  rotate  the  orbitals  as  the  calculations  require,  regardless  of  their  lack  of  optimum  property  for  a  given 
state.  Furthermore,  this  assists  in  the  ‘relative’  energies  of  the  multiplets  being  improved.  To  the  contrary,  using  different 
calculations  to  define  a  reference  complicate  the  interpretation  of  the  multiplet  separations,  since  different  orbitals  are  used  for 
different  reverence  functions. 


The  kind  of  results  we  obtain  are  shown  below.  The  difference  in  using  ‘fractional  occupation  number’  (FON)  orbitals,  which  are 
the  template  orbitals  for  all  these  systems,  differ  only  by  the  quantity  in  ()  from  results  using  the  most  optimum  set. 

FON 

CCSD  CCSD(T)  ^CCSD(T)  CCSDT-3  EXP 

Fe  5D(d6s2)  0.00  0.00  0.00  0.00  0.00 

5F(d7s1)  0.972 
(0.131)  1.068 
(-0.075)1.019 
(-0.044)0.949 
(0.005)  0.87 


Fe+6D(d6s1)  7.673 
(0.018)  7.846 
(-0.006)7.841 
(-0.013)7.839 
(-0.001)7.90 
4F(d7)  8.103 
(0.071)  8.365 
(-0.156)8.336 
(-0.144)8.194 
(-0.001)8.15 


Co  4F(d7s2)  0.00  0.00  0.00  0.00  0.00 

4F(d8s1)  0.480 
(0.201)  0.688 
(-0.154)0.616 
(-0.111)0.493 
(-0.014)0.42 


Co+  5F(d7s1)  8.066 

(0.01)  8.316 
(-0.087)8.325 
(-0.107)8.230 
(-0.002)8.28 
3F(d8)  8.055 
(-0.158)8.210 
(-0.305)8.167 
(-0.283)8.023 
(-0.138)7.85 


Ni  3D(d9s1)  -0.319 
(0.2751)  0.110 

(-0.278)0.001 
(-0.21 1)-0. 076 
(-0.059)-0.03 

3F  (d8s2)  0.00  0.00  0.00  0.00  0.00 

Ni+4F(sd8)  8.143 
(-0.006)8.349 
(0.001)  8.346 
(-0.009)8.460 
(-0.009)8.67 

2D(d9)  7.149 
(0.162)  7.504 
(-0.151)7.437 
(-0.115)7.455 
(-0.028)7.59 


Cu  2S  (dIOsI)  -1.554 
(0.346)  -0.966 
(-0.41 7)-1 . 1 1 2 
(-0.322)-1 .355 
(-0.104)-1.49 

2D(d9s2)  0.000  0.000  0.000  0.000  0.000 

2P(d10p1)  2.277 
(0.261)  2.711 
(-0.274)2.625 
(-0.227)2.315 
(0.00)  2.295 


Cu  + 

3D(d9s1)  8.828 
(-0.013)8.963 
(0.011)  8.967 
(0.00)  8.975 
(0.001)  9.04 

IS  (dIO)  6.036 
(0.213)  6.476 


(-0.239)6.392 
(-0.19)  6.255 
(-0.049)6.23 


Long  chain  hydrocarbons  and  their  folding. 

We  studied  the  issue  of  what  length  a  long-chain  hydrocarbon  can  turn  on  itself  to  form  a  hairpin  [3].  In  this  case,  we  studied 
alkanes  from  C8H18  to  C18H38.  The  accuracy  required  is  extreme  as  each  possible  gauch  rotation  from  the  standard  all  trans 
form  requires  only  2kj/mol.  The  calculations  were  only  made  possible  by  our  development  of  the  massively  parallel  ACES  III 
program  and  our  earlier  introduction  of  frozen  natural  orbitals  (FNO’s)  that  allows  us  to  eliminate  about  40%  of  the  virtual  space 
in  any  large  scale  correlated  calculation.  The  figures  attached  summarize  the  work.  All  geometries  were  optimized  at  the  CCSD 
level.  As  seen  from  the  figures  there  is  no  loss  in  accuracy  due  to  the  FNO’s. 

Fig  2 
Fig  3 
Fig  4 

Studies  of  nitramines  and  their  decomposition  pathways 
Fig  5 

Fig  6 
Fig  7 

The  mechanism  of  RDX  decomposition,  in  any  phase,  is  uncertain.  The  three  above,  NN  hemolysis,  HONO  elimination,  and  the 
so-called  ‘triple’whammy’  decomposition,  are  of  particular  interest.  Computational  prediction  of  the  gas-phase  mechanism  has 
been  hindered  by  methods  that  are  insufficiently  accurate  to  compare  energetically  close  mechanisms.  We  are  in  a  position  to 
describe  the  barrier  to  RDX  decomposition  in  the  gas  phase  with  unprecedented  accuracy  using  coupled  cluster  theory,  the 
computational  method  against  which  other  methods,  like  DFT,  are  compared  against  for  reliability. 

We  have  recently  reported  the  electronic  and  free  energy  barriers  to  RDX  decomposition  using  CCSD(T)/CBS  (complete  basis 
set)  electronic  energies  with  MBPT(2)/cc-pVTZ  structures  and  partition  functions  using  Eyring  transition  state  theory.  We 
investigate  all  major  considered  mechanisms  for  decomposition  and  subsequent  decomposition  pathways  associated  therein, 
including  NN  homolysis,  HONO  elimination,  triple  whammy,  and  NONO  isomerization.  We  find  that  the  dominant  mechanism  is 
HONO  elimination;  there  is  not  a  competition  between  HONO  elimination  and  NN  homolysis,  as  DFT  predictions  previously 
predicted.  We  further  show  the  adjustability  of  the  HONO  barrier  to  chemical  substitution,  providing  a  proof-of-principle  of  how 
one  may  adjust  RDX,  chemically,  to  change  its  shock-sensitivity  value. 

The  decomposition  mechanism  of  RDX  has  been  studied  for  about  70  years  without  a  decisive  specific  conclusion.  There  are 
many  challenges  to  establishing  the  mechanism.  It  is  currently  not  possible  to  reconcile  all  the  conclusions  from  different 
experiments,  as  they  are  mutually  exclusive.  Some  of  the  challenges  to  analysis  include: 

a)  RDX  has  closely  related  molecular  cousins,  HMX  and  CL-20.  To  what  extent  do  the  experimental  results  of  these  very 
similar  compounds  apply  to  RDX? 

b)  In  which  phase  does  the  reaction  begin?  Evidence  exists  for  solid,  liquid,  and  gas. 

c)  Independent  of  the  phase  in  which  the  reaction  begins,  in  which  phase  does  the  chemical  reaction  most  proceed? 

d)  Does  the  evidence  gathered  by  a  specific  experimental  mean  relate  to  the  “natural”  RDX  decomposition?  It  is  possible  that 
some  of  the  evidence  is  not  in  contradiction,  but  that  mass  spectrometry  vs.  UV-VIS  excitation  vs.  dropping  an  anvil  on  RDX  vs. 
heating  all  inherently  induce  different  mechanisms. 

e)  Which  small-molecule  gaseous  products  are  observed  first?  Studies  focus  on  C02,  N02,  NO,  N20,  H20,  HCN,  and  many 
others  as  clues  to  the  mechanism. 

f)  Do  multiple  mechanisms  occur  simultaneously? 

g)  Is  the  mechanism  of  initiation  the  dominant  mechanism  of  propagation? 

h)  What  computational  methods  have  the  accuracy  to  estimate  energy  barriers  of  transition  states  accurately? 

i)  Independent  of  speculating  on  quasi  non-observables  (the  “mechanism”  itself  is  not  an  observable),  how  does  one 
reconcile  the  extreme  disagreement  in  clear  observables,  like  activation  energy? 

We  are  in  a  unique  position  to  provide  clarity  to  a  complex  problem.  Computational  modelling  provides  unambiguous  results 
in  terms  of  meaning;  we  do  not  need  to  divine  complex  phenomena  from  isolated  points  in  time  of  one  physical  property  alone. 
There  is  no  speculation  to  be  had  about  what  is  being  said.  We  are,  however,  restricted  in  the  realism  of  our  modelling. 
Modelling  liquid  or  solid  systems  requires  simulation  of  many  particles;  modelling  of  the  gas-phase,  however,  is  very 
straightforward.  It  is  our  intention  to  solve  the  gas-phase  decomposition  of  RDX. 

This  has  three  purposes: 

1)  In  this  way,  experimentalists  may  gauge  which  techniques  reproduce  the  gas-phase  results,  and  then  apply  those 
techniques  to  the  much  more  complex  liquid  and  solid  phases. 

2)  To  a  first-order  approximation,  condensed-phase  mechanisms  are  well  described  by  gas-phase  kinetics.  There  are  some 
notable  exceptions  in  which  the  potential  energy  surface  is  qualitatively  different  between  condensed  and  gas-phase,  but  at 


worse,  the  gas-phase  mechanism  provides  significant  guidance  (especially  for  a  molecular  crystal,  whose  intermolecular  forces 
are  weak). 

3)  In  establishing  the  mechanism,  we  have  insight  into  how  to  may  adjust  the  shock-sensitivity  of  RDX  and  its  HMX  and  CL20 
cousins. 
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The  ACES  III  program  system  being  developed  un  this  grant  is  in  wide  use  at  the  DOD  HPCC.  Active  users  at  ARL  include  de 
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SR-CC,  V=exp(T)<P 

BASIC  MR  PROBLEM:  Occurs  when  orbs 
like  I  and  A  become  equally  important. 
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C0-  90%  Defines  SRCC  problem. 
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When  the  C’s  are  too  large  to  achieve  this  value  via  T 
in  an  applicable  SRCC  method,  we  encounter  a  MR-CC 
problem.  We  need  a  way  to  obtain  the  C's  without  the  SR  bias. 


J  _T_i_ 


Fig  2 


Illustration  of  an  accurate  energy  problem:  At  what  length 
does  dispersion  cause  normal  alkanes  to  assume  a  hairpin  form? 


The  Last  Globally  Stable  Extended  Alkane*' 

Nils  O  B  Lttschwager,  Tobias  N,  Wassermann, 
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Angew  Chem,  Int.  2012. 

Exp  points  to  Cl 8. 

Each  gauch  rotation  makes  -  2kJ/mol  difference 
Jason  Byrd,  RJB  J  Montgomery  ,  JPC  2014 
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n-alkanes,  CpH2p*2 

electrons  basis  functions 

C8 

66 

492 

CIO 

82 

608 

C12 

98 

724 

C14 

114 

840 

C16 

130 

956 

Cl  8 

146 

1072 

C12  29  hrs.  CCSD(T)/vtz  489  processors 

6  hrs.  FNO(CCSD(T))/vtz  240  processors 

C14  70  hrs.  CCSD(T)/vtz  832  processors 

13  hrs.  FNO(CCSD(T))/vtz  480  processors 
C18  51  hrs.  FNO(CCSD(T))/vtz  960  processors 
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At  what  alkane  length  Is  the  hairpin  form  preferred  over  the  linear  one? 


Fig  5 
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